C*************************************************************** C Test program for subroutine LSEM3 by Stephen Kirkup * C*************************************************************** C C Version 2 (July 2015) (Version 1 2001) Changes in LIBEM3. C School Engineering, University of Central Lancashire C www.uclan.ac.uk C smkirkup@uclan.ac.uk C http://www.researchgate.net/profile/Stephen_Kirkup C C This open source code can be found at C www.boundary-element-method.com/fortran/LSEM3_T.FOR C C Issued under the GNU General Public License 2007, see gpl.txt C C Part of the the author's open source BEM packages. C All codes and manuals can be downloaded from C www.boundary-element-method.com C C*************************************************************** C This program is a test for the subroutine LSEM3. The program computes C the solution to a Laplace problem exterior to two square plates by C the shell panel method. One plate is assigned the potential one C and the other is at zero potential. The problem may be thought of as C a pair of plates of a capacitor with one plate at zero volts, the C other at 1 volt. C C Background C ---------- C C We wish to solve the Laplace equation C C __ 2 C \/ {\phi} = 0 C C C For the exterior problem, the domain lies exterior to the shells C or plates. The shell condition may be Dirichlet, Robin or C Neumann. It is assumed to have the following general form C C a(q) {\delta}(q) + b(q) {\nu}(q) = f(q) C C A(q) {\Phi}(q) + B(q) V(q) = F(q) C C where {\delta}(q)={\phi+}(q) - {\phi-}(q), C {\nu}(q) = {\nu+}(q) - {\nu-}(q) C {\Phi}(q) = ({\phi+}(q)+{\phi-}(q))/2 C V(q) = ({v+}(q) + {v-}(q))/2 C and {\phi}(q) is the velocity potential at the point q on S, v(q) C is the derivative of {\phi} with respect to the outward normal to S C at q and {\alpha}, {\beta} and f are real-valued functions defined C on S. C C Subroutine LSEM3 accepts a description of the shells and the position C of the exterior points where the solution ({\phi}) C is sought, the shell condition and returns the solution ({\delta}, C {\Phi}, {\nu} and V) on {\Pi} and the value of {\phi} at the exterior C points. C C The test problems C ----------------- C C In this test the domain is the exterior to two square plates of side C 0.1m and 0.01m apart. The solution to the problem with a Dirichlet C shell condition (a=1, b=0, f=0) on both plates and (A=1, B=0, C F=0) on the lower plate and (A=1, B=0, F=1) on the upper plate. C C Each plate is described by a set of 32 planar triangular panels C of equal size. The shell solution points are the centres of the C panels. The solution is sought at the exterior points C (0.05,0.05,z), the axis through the plates with z=0.001,0.003, C 0.005, 0.007, 0.009. C---------------------------------------------------------------------- C The PARAMETER statement C ----------------------- C There are four components in the PARAMETER statement. C integer MAXNH : The limit on the number of shell panels. C integer MAXNV : The limit on the number of vertices. C integer MAXNPE : The limit on the number of exterior points. C External modules related to the package C --------------------------------------- C subroutine LEBEM3: Subroutine for solving the exterior Laplace C equation (file LEBEM3.FOR contains the subroutine LEBEM3) C subroutine H3LC: Returns the individual discrete Laplace integral C operators. (file H3LC.FOR contains H3LC and subordinate routines) C subroutine GLS: Solves a general linear system of equations. C (file GLS.FOR contains CGSL and subordinate routines) C file GEOM3D.FOR contains the set of relevant geometric subroutines C The program PROGRAM LSEM3T IMPLICIT NONE C VARIABLE DECLARATION C -------------------- C PARAMETERs for storing the limits on the dimension of arrays C Limit on the number of panels INTEGER MAXNH PARAMETER (MAXNH=10) C Limit on the number of vertices (equal to the number of panels) INTEGER MAXNV PARAMETER (MAXNV=12) C Limit on the number of points exterior to the plate, where C potential properties are sought INTEGER MAXNPE PARAMETER (MAXNPE=5) C Geometrical description of the plate(s) C Number of panels and counter INTEGER NH,IPI C Number of vetices and counter INTEGER NV,IV C Index of nodal coordinate for defining boundaries (standard unit is C metres) REAL*8 VERTEX(MAXNV,3) C The three nodes that define each panel on the boundaries INTEGER HELV(MAXNH,3) C The points exterior to the plate(s) where the potential C properties are sought and the directional vectors at those points. C [Only necessary if an exterior solution is sought.] C Number of exterior points and counter INTEGER NPE C Coordinates of the exterior points REAL*8 PEXT(MAXNPE,3) C Shell Conditions C a(p)+b(p)=f(p) C function a at the central points REAL*8 HA(MAXNH) C function b at the central points REAL*8 HB(MAXNH) C function f at the central points REAL*8 HF(MAXNH) C A(p)+B(p)=F(p) C function A at the central points REAL*8 HAA(MAXNH) C function B at the central points REAL*8 HBB(MAXNH) C function F at the central points REAL*8 HFF(MAXNH) C Incident field C The incident potential {\phi} at the central points REAL*8 HIPHI(MAXNH) C The derivative of the incident potential {\phi} at the C central points REAL*8 HIVEL(MAXNH) C The incident potential {\phi} at the exterior points REAL*8 PEIPHI(MAXNPE) C Validation and control parameters LOGICAL LSOL LOGICAL LVALID REAL*8 EGEOM C Solution REAL*8 PHIDIF(MAXNH) REAL*8 PHIAV(MAXNH) REAL*8 VELDIF(MAXNH) REAL*8 VELAV(MAXNH) REAL*8 PEPHI(MAXNPE) C Working space REAL*8 AMAT(2*MAXNH,2*MAXNH) REAL*8 BMAT(2*MAXNH,2*MAXNH) REAL*8 L_EH(MAXNPE,MAXNH) REAL*8 M_EH(MAXNPE,MAXNH) C Further information from system solution GLS INTEGER*4 PERM(2*MAXNH) LOGICAL XORY(2*MAXNH) REAL*8 C(2*MAXNH) C Work Space REAL*8 WKSPC1(2*MAXNH) REAL*8 WKSPC2(2*MAXNH) REAL*8 WKSPC3(2*MAXNH) C Counter through the x,y,z coordinates INTEGER ICOORD C Counter through the exterior points INTEGER IPE C Distance between plates REAL*8 GAP C INITIALISATION C -------------- C Describe the nodes and panels that make up the plate C :The two 1x1 plates are each divided into 32 panels, C : (NH=64). VERTEX and HELV are defined anti-clockwise around C : the plate when viewed from above so that the normal to the C : plate is assumed to be upward C :Set up nodes C : Set NH, the number of panels NH=8 C : Set NV, the number of vertices (equal to the number of panels) NV=10 C : Set coordinates of the nodes C : Set up VERTEX, the coordinates of the vertices of the panels DATA ((VERTEX(IV,ICOORD),ICOORD=1,3),IV=1,10) C Plate 0 (z=0) * / 0.000D0, 0.000D0, 0.000D0, * 1.000D0, 0.000D0, 0.000D0, * 1.000D0, 1.000D0, 0.000D0, * 0.000D0, 1.000D0, 0.000D0, * 0.500D0, 0.500D0, 0.000D0, C Plate 1 (z=1) * 0.000D0, 0.000D0, 0.100D0, * 1.000D0, 0.000D0, 0.100D0, * 1.000D0, 1.000D0, 0.100D0, * 0.000D0, 1.000D0, 0.100D0, * 0.500D0, 0.500D0, 0.100D0/ C : Set nodal indices that describe the panels of the shells. C : The indices refer to the nodes in VERTEX. The order of the C : nodes in HELV dictates that the normal is outward from the C : plate into the potential domain. DATA ((HELV(IPI,ICOORD),ICOORD=1,3),IPI=1,8) C Plate 0 (z=0) * / 1, 2, 5, 2, 3, 5, 3, 4, 5, 4, 1, 5, C Plate 1 (z=0.5) * 6, 7, 10, 7, 8, 10, 8, 9,10, 9, 6, 10/ C Set NPE NPE=5 C Set coordinates of chosen exterior points GAP=0.1D0 DO 300 IPE=1,NPE PEXT(IPE,1)=0.5D0 PEXT(IPE,2)=0.5D0 PEXT(IPE,3)=GAP*DFLOAT(2*IPE-1)/DFLOAT(2*NPE) 300 CONTINUE C Set the potential on the plates DO 310 IPI=1,NH HA(IPI)=1.0D0 HB(IPI)=0.0D0 HF(IPI)=0.0D0 HAA(IPI)=1.0D0 HBB(IPI)=0.0D0 C Plate 1 IF (IPI.LE.4) HFF(IPI)=0.0D0 C Plate 2 IF (IPI.GT.4) HFF(IPI)=1.0D0 310 CONTINUE DO 320 IPI=1,NH HIPHI(IPI)=0.0D0 HIVEL(IPI)=0.0D0 320 CONTINUE DO 330 IPE=1,NPE PEIPHI(IPE)=0.0D0 330 CONTINUE C :Switch for particular solution LSOL=.TRUE. C :Switch on the validation of LSEM3 LVALID=.TRUE. C :Set EGEOM EGEOM=1.0D-6 CALL LSEM3(MAXNV,NV,VERTEX,MAXNH,NH,HELV, * MAXNPE,NPE,PEXT, * HA,HB,HF,HAA,HBB,HFF, * HIPHI,HIVEL,PEIPHI, * LSOL,LVALID,EGEOM, * PHIDIF,PHIAV,VELDIF,VELAV,PEPHI, * AMAT,BMAT,L_EH,M_EH, * PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3) 999 FORMAT(3F14.4) OPEN(10,FILE='LSEM3.OUT') WRITE(10,*) 'LSEM3: Test Problem Results' WRITE(10,*) ' z computed pot exact/app pot' DO 600 IPI=1,NPE WRITE(10,999) PEXT(IPI,3),PEPHI(IPI),PEXT(IPI,3)*10 600 CONTINUE END